Известны ежемесячные продажи австралийского вина в тысячах литров с января 1980 по июль 1995, необходимо построить прогноз на следующие два года.

Попробуем поделить на число дней в месяце: Ряд не стал более регулярным, так что вернёмся к исходным данным.

STL-декомпозиция ряда:

Оптимальное преобразование Бокса-Кокса и результат его применения:

В данном случае преобразование имеет смысл использовать, так как оно хорошо стабилизирует дисперсию. Попробуем округлить параметр и взять \(\lambda=-1\): Результат практически такой же. Далее будем использовать \(\lambda=-1\).

Прогноз ETS

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = tSeries, lambda = LambdaOpt) 
## 
##   Box-Cox transformation: lambda= -1 
## 
##   Smoothing parameters:
##     alpha = 0.0707 
##     beta  = 0.0043 
##     gamma = 1e-04 
##     phi   = 0.9782 
## 
##   Initial states:
##     l = 0.9999 
##     b = 0 
##     s=0 0 0 0 0 0
##            0 0 0 0 0 0
## 
##   sigma:  0
## 
##       AIC      AICc       BIC 
## -3453.752 -3449.368 -3396.786

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set   51.61557 2189.653 1642.288 -0.3744397 6.478496 0.8362138
## Test set     1232.80660 2127.679 1766.955  3.6014104 6.978830 0.8996913
##                     ACF1 Theil's U
## Training set -0.08120999        NA
## Test set     -0.35511396 0.3246609

Остатки:

Достигаемые уровни значимости критерия Льюнга-Бокса для них:

Остатки не коррелированы

Q-Q plot и гистограмма для остатков:

Распределение имеет выброс слева.

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка отвергается 1.900071910^{-4}
Несмещённость Уилкоксона не отвергается 0.6891151
Стационарность KPSS не отвергается 0.1

ARIMA

Автоматический подбор модели

Применим функцию auto.arima:

## Series: tSeries 
## ARIMA(1,1,2)(0,1,2)[12] 
## Box Cox transformation: lambda= -1 
## 
## Coefficients:
##           ar1      ma1      ma2     sma1     sma2
##       -0.5897  -0.2604  -0.5390  -0.4770  -0.1787
## s.e.   0.6409   0.6149   0.5321   0.0918   0.0866
## 
## sigma^2 estimated as 8.281e-08:  log likelihood=1780.08
## AIC=-3548.15   AICc=-3547.61   BIC=-3529.63

Предлагается модель ARIMA(1,1,2)(0,1,2)\(_{12}\). Посмотрим на её остатки:

Отрежем первые 13 отсчётов и продолжим анализ:

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка отвергается 6.367987210^{-5}
Несмещённость Уилкоксона отвергается 0.0460378
Стационарность KPSS отвергается 0.0327717

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:

##                     ME     RMSE      MAE        MPE      MAPE      MASE
## Training set 543.62191 4285.153 2647.961  2.7109642 11.575578 1.3482784
## Test set     -56.43371 1961.410 1517.384 -0.9025861  6.206568 0.7726154
##                    ACF1 Theil's U
## Training set  0.6188588        NA
## Test set     -0.3442693 0.2790998

Ручной подбор модели

Исходный ряд нестационарен (p<0.01, критерий KPSS); сделаем сезонное дифференцирование: Ряд всё ещё нестационарен (p<0.01, критерий KPSS). Проведём ещё одно дифференцирование: Для полученного ряда гипотеза стационарности не отвергается (p>0.1)

Посмотрим на ACF и PACF полученного продифференцированного ряда:

На ACF значимы лаги 1, 12 и 24, на PACF — 1-5. Поищем с помощью auto.arima оптимальную модель полным перебором (stepwise=FALSE) с ограничениями d=1, D=1, max.p=5, max.q=4, max.P=0, max.Q=2, max.order=10:

## Series: tSeries 
## ARIMA(1,1,1)(0,1,2)[12] 
## Box Cox transformation: lambda= -1 
## 
## Coefficients:
##           ar1      ma1     sma1    sma2
##       -0.0201  -0.8693  -0.4843  -0.177
## s.e.   0.0838   0.0415   0.0923   0.087
## 
## sigma^2 estimated as 8.229e-08:  log likelihood=1780.02
## AIC=-3550.05   AICc=-3549.66   BIC=-3534.61

Была найдена более оптимальная по AICc модель — ARIMA(1,1,1)(0,1,2)\(_{12}\). Посмотрим на её остатки: Отрежем начало ряда остатков и проанализируем их:

Достигаемые уровни значимости критерия Льюнга-Бокса для остатков:

Q-Q plot и гистограмма:

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка отвергается 4.503317210^{-5}
Несмещённость Уилкоксона не отвергается 0.055308
Стационарность KPSS отвергается 0.0382442

Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:

##                      ME     RMSE      MAE        MPE      MAPE      MASE
## Training set 535.674930 4312.632 2699.617  2.6262660 11.786210 1.3745806
## Test set       4.690032 1936.939 1490.400 -0.7019034  6.133156 0.7588758
##                    ACF1 Theil's U
## Training set  0.6172789        NA
## Test set     -0.3603302 0.2778012

Сравним остатки двух версий аримы, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены: 

## 
##  Diebold-Mariano Test
## 
## data:  resres.auto
## DM = 0.21196, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.8324
## alternative hypothesis: two.sided

Критерий Диболда-Мариано не обнаруживает значимого различия между качеством прогнозов.

В целом подобранная вручную модель проще, её остатки лучше, а ошибка на тесте меньше, так что остановимся на модели, подобранной вручную.

Итоговое сравнение

Сравним остатки лучших моделей ARIMA и ETS, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены:

## 
##  Diebold-Mariano Test
## 
## data:  res.autores.ets
## DM = 1.4329, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.1538
## alternative hypothesis: two.sided

Критерий Диболда-Мариано не обнаруживает значимого различия между качеством прогнозов, но у ARIMA ошибка на тесте и AIC меньше, так что выберем её.

Финальный прогноз

Поскольку остатки ненормальны, предсказательные интервалы строятся бутстрепом

## Warning in InvBoxCox(pred$pred, lambda, biasadj, var(residuals(object), :
## biasadj information not found, defaulting to FALSE.
##          Point Forecast    Lo 80 Hi 80    Lo 95 Hi 95
## Aug 1994       29266.72 2488.820    NA 1676.706    NA
## Sep 1994       24226.65 2432.219    NA 1647.597    NA
## Oct 1994       27958.44 2447.141    NA 1650.091    NA
## Nov 1994       32727.79 2460.692    NA 1651.952    NA
## Dec 1994       37983.85 2468.522    NA 1651.224    NA
## Jan 1995       15250.22 2235.830    NA 1540.085    NA
## Feb 1995       21981.02 2325.168    NA 1578.128    NA
## Mar 1995       24016.60 2330.496    NA 1576.790    NA
## Apr 1995       25868.30 2331.197    NA 1573.366    NA
## May 1995       24155.14 2301.301    NA 1556.054    NA
## Jun 1995       25890.01 2301.108    NA 1552.372    NA
## Jul 1995       30204.94 2315.639    NA 1555.389    NA
## Aug 1995       28523.51 2004.960    NA 1343.665    NA
## Sep 1995       24887.50 1965.200    NA 1321.083    NA
## Oct 1995       27791.14 1960.512    NA 1313.995    NA
## Nov 1995       32830.91 1961.203    NA 1309.437    NA
## Dec 1995       38462.01 1958.033    NA 1303.254    NA
## Jan 1996       15994.05 1810.062    NA 1231.788    NA
## Feb 1996       21658.53 1847.446    NA 1244.730    NA
## Mar 1996       24290.08 1847.047    NA 1240.366    NA
## Apr 1996       25943.91 1838.750    NA 1232.530    NA
## May 1996       24453.66 1814.179    NA 1217.492    NA
## Jun 1996       25509.52 1803.462    NA 1208.801    NA
## Jul 1996       30508.44 1808.362    NA 1207.192    NA